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Anisotropic exciton Stark shift in black phosphorus 
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We calculate the excitonic spectrum of few-layer black phosphorus by direct diagonalization of the 
effective mass Hamiltonian in the presence of an applied in-plane electric field. The strong attractive 
interaction between electrons and holes in this system allows one to investigate the Stark effect up to 
very high ionizing fields, including also the excited states. Our results show that the band anisotropy 
in black phosphorus becomes evident in the direction dependent field induced polarizability of the 
exciton. 
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Excitons in semiconductors have been subject of inves¬ 
tigation for many years. Qi Such interacting electron- 
hole pair mimics a hydrogen atom, with the hole play¬ 
ing the role of the nucleus. i It is well known that 
the hydrogen atom in the presence of an applied elec¬ 
tric field exhibits degeneracy breaking and a quadratic 
energy shift due to the so called Stark effect. The exper¬ 
imental observation of such an effect for excitons in bulk 
semiconductors is however hindered by the low electron- 
hole binding energy - for GaAs, for example, this energy 
is around 4.8 meV, Q allowing only for very low ioniz¬ 
ing electric fields. In such low fields, the exciton binding 
energy is not significantly modified, albeit the exciton 
peak broadens due to the decrease of the exciton life¬ 
time. This is a consequence of exciton ionization, which 
makes the experimental observation of the exciton Stark 
effect much more challenging. This motivates the study 
of the Stark effect in quantum wells, which, through the 
so-called quantum confined Stark effect, can circumvent 
electron-hole dissociation, allowing the use of higher elec¬ 
tric fields. 0 Alternatively, the excitonic Stark effect has 
been also theoretically investigated in carbon nanotubes, 
[^, 0 where the electron-hole binding energies, depend¬ 
ing on the nanotube configuration, may reach quite large 
values, [H allowing for high ionizing fields. Nevertheless, 
setting up an experiment to detect this effect in carbon 
nanotubes is a difficult task, which has been achieved 
only very recently. @ Strong exciton binding energies are 
also observed in conjugated polymer chains, 0 where 
exciton Stark shifts for hig h electric fields have been in¬ 
vestigated as well. [ll|, ll2| 

In recent sudies on single or few-layer semiconductors, 
such as transition metal dichalcogenides and black phos¬ 
phorus, exciton binding e nerg ies are found to be on the 
order of hundreds of meV, |l3l4l6j which brings the possi¬ 
bility of experimentally observing the Stark effect of their 
excitons. In this context, the case of black phosphorus, 
a layered material that has recently been fabri¬ 
cated in few-layer form 0 and has a strong potential for 


technological applications, lH, [22l - [^ is of special inter¬ 
est. Its effective mass anisotropy |27H29] leads to an exci¬ 
ton wave function with distinct distributions in different 
in-plane directions, so that the exciton Stark shift behav¬ 
ior, namely, its electric field induced dipole moments and 
polarizability, is expected to depend on the direction of 
the applied in-plane electric field, as we will demonstrate 
by our numerical calculations in what follows. 

We consider the exciton Hamiltonian within the effec¬ 
tive mass approximation 


Hexc — Hq + Hy E(|re 


( 1 ) 


Due to the high anisotropy of black phosphorus energy 
bands, it is convenient to write the single particle Hamil¬ 
tonian in cartesian coordinates {xe^ye^Xh^Vh) as 


Hi = - 




2mix dx‘1 2miy dyf 


+ QiF • fl 


( 2 ) 


where i = e{h) represents the electron (hole), F is the 
in-plane applied electric field, rriix^y) is the effective mass 
of the carrier i in the x(^)-direction, and qi is its charge. 
It is convenient to rewrite the Hamiltonian in units of 
the Rydberg energy Ry and Bohr radius ao. We then use 
relative (x, y) and center-of-mass (R) coordinates in each 
direction to simplify the exciton Hamiltonian as 


1 

Ex dx‘^ 


Y^+V{^T^^TT)+eF-r, (3) 
Ty uy 

where Rx^y) is the reduced effective mass in the x{y)- 
direction and the center-of-mass contribution is removed, 
since the potential does not depend on these coordinates. 
In order to take proper account of the large dielectric 
contrast between the vacuum on top of the sample and 
the substrate below it, the interaction potential E(r), 
for r = is assumed to be of the Keldysh type 


2^, [ 3 ^ , which in dimensionless form reads 


V{r) = - 


27r 


(ei + e2)po 


Ho{—)-Yo(- 
PoJ \Po/i 


(4) 
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where Hq and Iq are Struve and Neumann functions, 
respectively, ei( 2 ) is the dielectric constant of the vac¬ 
uum (substrate) surrounding the phosphorene layer, and 
Po = De /(ei + 62 ) is the screening factor, where D is the 
effective width of the layer and e is its dielectric constant, 
which is isotropic for our problem of interest, as discussed 
in Ref. l3l|. For ni layers of phosphorene on a Si 02 sub¬ 
strate, this factor was found to be po = n/10.79A. 0 
To find the excitonic eigenstates, one can use varia¬ 
tional functions as proposed in Refs. HI, [HI- However, 
this kind of approach normally allows one to find only the 
ground state energy, or just a few excited states. Hence, 
in order to obtain several exciton states, we rather nu¬ 
merically diagonalize H within a finite difference scheme, 
which must be performed with a variable mesh, in or¬ 
der to take better account of the singularity of V (r) as 
r ^ 0. Technical details of the finite difference approach 
in a variable mesh are provided in the supplementary 
material. [33 



FIG. 1: Energies of low-lying exciton states, considering a 
black phosphorus sample with nz = 1 (black circles), 2 (red 
triangles), 3 (blue squares) and 4 (green stars) phosphorene 
layers. Fitting functions for each case are shown as curves. In¬ 
set: exciton energy states for WS 2 (black squares) and for an 
electron-hole pair interacting by the Coulomb potential (red 
circles) with an effective relative permittivity Cr, along with 
experimentally obtained values (blue stars, courtesy from the 
authors of Ref. 0 ), for comparison. Degenerate states are 
connected by horizontal lines, to help visualization. 

The obtained exciton energies are shown in Fig. [T]for 
four values of thickness (number of phosphorene layers) of 
the black phosphorus sample, calculated using the same 
effective masses as in Ref. [l3|. For the sake of com¬ 
parison, the inset shows (i) the results for WS 2 , where 
conduction and valence bands are known to be isotropic 
with /i = /ix = liy = 0.16; [l5[ (ii) results for the case 
where the electron and hole interact via the Coulomb 
potential, i.e. for the 2D hydrogen atom, where the per¬ 
mittivity e ~ 5.23eo is adjusted as to fit the ground state 


energy of WS 2 , by means of the analytical expression 

= ARyiijE^^ where Eq = 0.318 eV is the ground state 
exciton binding energy; and (iii) experimental results. 

0 

Let us first discuss the results in the inset. Notice 
that the energies obtained from the Coulomb potential 
are clearly underestimated as compared to those calcu¬ 
lated with the Keldysh potential. This is a consequence 
of the fact that the former decays faster than the latter as 
r ^ 0. Moreover, results from the Coulomb potential are 
analytically found to follow En = —Ryl{n — ll2)‘^^ where 
the states are (2n+l)-fold degenerate. In fact, this is eas¬ 
ily verified with the numerical method proposed in this 
work. We point out that binding energies calculated with 
the actual permittivity of WS 2 , e = 13eo, would be even 
more underestimated, since the effective Rydberg energy 
is inversely proportional to e^. On the other hand, taking 
advantage of the fact that the potential in Eq. (|4j) has 
a 1/r tail for r ^ oc, one can adjust e 1.8eo as to fit 
the numerically obtained higher exciton energy states - 
in this case, due to such smaller effective e, the Coulomb 
interaction model provides highly overestimated binding 
energies for the lower states. A two-fold degeneracy is 
still expected due to the circular symmetry of the system, 
so that positive or negative values of angular momentum 
I lead to the same energy. However, in the case of the 
Coulomb potential, additional degeneracies are also ob¬ 
served, due to the SO(3) symmetry of this potential in 
two dimensions, which is not the case for the Keldysh 
potential for WS 2 . Note that for black phosphorus, not 
even the two-fold degeneracies are observed, since the an¬ 
gular momentum is not a good quantum number in this 
case, due to the anisotropy of the kinetic energy. 

Inspired by the analytical expression for the eigenstates 
of the Coulomb potential and regarding the fact that 
the Keldysh potential in Eq. 0 has the Coulomb form 
as asymptotic behavior, we fit our numerical results by 
En = — 7 i/(n + 72 )^^, as shown by the curves in Fig. [H 
Fitting parameters are given in Table I. The ground state 
exciton energies for n/ =1, 2, 3 and 4 phosphorene layers 
are found to be -396, -261, -200 and -163 meV, respec¬ 
tively. This is in good agreement with (namely, less than 
10% lower than) the results obtained by the variational 
approach. 0 The three lowest energy s states of WS 2 , 
are found to be -318.80, -151.72, and -93.56 meV, also in 
good agreement with recent experimental results for this 
material, which are shown as open blue stars in the inset 
of Fig. [H Q which validates our method. 

By using Fermi’s golden rule, one can calculate the 
probability for a photon induced valence band-to-exciton 
transition to occur, which is demonstrated to be propor¬ 
tional to |(c|e • ^'r’)p|' 0 n(O 5 0 )p 5 where the first factor is 
the dipole matrix element between conduction and va¬ 
lence band states. The optically active transition in the 
case of light polarized in the x—direction is between the 
band states labeled as F^ and F^ in Ref. l28|, which is 
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TABLE I: Fitting parameters for the exciton spectra in Fig. [U 
using the Rydberg-like expression En = 7 i/(n+ 72 )"’'^, consid¬ 
ering different number of phosphorene layers. The parameter 
7 i is in units of meV, whereas 72 and 73 are dimensionless. 


ni 

1 

2 

3 

4 

71 

317.38 

223.92 

183.51 

170.45 

72 

-0.24 

-0.19 

-0.119 

-0.0586 

73 

-0.80 

-0.72 

-0.699 

-0.733 


the one we deal with here. In fact, such a polarization 
in the lighter effective mass direction has been recently 
verified experimentally. [HI With this factor set, we just 
need now to investigate the squared modulus of the exci¬ 
ton envelope function at r = 0, i.e. ^ as being 

the oscillator strength. Indeed, it is reasonable that elec¬ 
tron hole recombination processes are more likely when 
electron and hole are in the same place in real space, i.e., 
if |'0^(O,O)p is non-zero. Therefore, for isotropic mate¬ 
rials such as WS 2 , the optically active exciton states are 
only those with 5 -like wavefunctions. This means that 
the exciton Rydberg series observed in the experiments 
performed in Ref. [l^ are not the complete exciton en¬ 
ergy series: I 7 ^ 0 states are missing in the spectrum, 
and only the states labelled as I 5 , 25 and 35 in the inset 
of Fi g. [J are observed experimentally in e.g. Fig. 3 of 
Ref. [ 1 ^, with very good quantitative agreement. The 
p states would appear in two-photon experiments, yet 
with a different energy series, in contrast with the case 
of Coulomb-like electron-hole interaction, where for each 
p state there is a degenerate 5 state. In fact, the p states 
in between I 5 and 25 states observed in our results are 
consistent with recent experimental observations, In 
the case of black phosphorus, the angular momentum is 
no longer a good quantum number, since the kinetic en¬ 
ergy operator is not circularly symmetric, but we observe 
that n = 1 {El = -396.37 meV) and n = 3 {E^ = -142.78 
meV) states in the spectrum of Fig. [T] are optically ac¬ 
tive, i.e. |'0i(O,O)p 7 ^ 0 and |'03(O,O)p 7 ^ 0, due to 5 -like 
orbital nature of their wave functions (details of the wave 
functions for these and other exciton states are shown in 
the supplementary material). Thus, we will refer only 
to these states in the discussion about the exciton Stark 
effect from here onwards. 

Due to the anisotropic band structure of black phos¬ 
phorus, some physical properties of this material are ex¬ 
pected to be direction dependent. Exciton response to 
applied electric fields, for instance, should be stronger in 
the direction where the reduced effective mass is lower. 
This is indeed observed in Fig. O which shows, as sym¬ 
bols, the numerically obtained Stark shift for an electric 
field applied in x (left panels) and ^-directions (right pan¬ 
els), considering n/ = 1 - 4 layers. The quadratic Stark 



FIG. 2 : Ground (n = 1 ) state exciton stark shift, consider¬ 
ing electric fields applied in x (left panels) and ^-direction 
(right panels), for systems with ni = 1 (black circles), 2 (red 
triangles), 3 (blue squares) and 4 (green stars) phosphorene 
layers. Expectation values of the electron-hole polarization 
as a function of the applied electric field intensity are shown 
in the bottom panels. Gurves in top and bottom panels are 
quadratic and linear fittings for the stark shift and polariza¬ 
tion, respectively. 


shift follows the expression 

AE = E{F) - E{E = 0)=pF + I3F^ (5) 

where F is the magnitude of the electric field, p is the 
intrinsic dipole moment, and /3 is the polarizability pa¬ 
rameter. Such a quadratic behavior is, to a good extent, 
observed in our numerical results, which are well fitted by 
the quadratic curves in the top panels of Fig. O where we 
assumed p = 0 and used /3 as a fitting parameter in Eq. 
([5]). The p = 0 assumption is justified by the fact that 
we are dealing only with non-degenerate states, which 
do not possess a permanent dipole moment, so that the 
linear Stark shift can be neglected. The quadratic Stark 
shift comes from the fact that the applied field is able to 
induce electron-hole polarization, which depends linearly 
on the strength of the applied field. Such linear behav¬ 
ior is indeed observed in the bottom panels of Fig. [2l 
where the numerically obtained polarizations (symbols), 
i.e. the expectation values of the electron-hole separation 
{x) and (y), are well fitted by linear functions (curves). 
As observed in Fig. [TJ the ground state exciton energy 
decreases as the number of layers increases. As a con¬ 
sequence, the exciton energy becomes more succeptible 
to electric field effects for higher number of layers and 
the Stark shift is more clear in these cases - for instance, 
considering n/ = 4, the exciton ground state exhibits a 
Stark shift of ~ 4.5 meV and a polarization of {x) ~ 25 
A for fields as high as F = 100 kV/cm applied in the x 
direction, whereas a much smaller AE ^ 1.5 meV and 
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TABLE II: Polarizability parameters and l3y for the n = 1 
and 3 states, in units of eA^/mV, for electric fields applied in 
the X and y direction, respectively. 


ni 

1 

2 

3 

4 

Px{y) 

0.37 (0.135) 

1.06 (0.42) 

2.05 (0.93) 

3.5 (1.5) 

Pl(y) 

10 (4.8) 

27 (13) 

48 (23) 

65 (43) 


{x) ^ 4.5 A is observed for n/ = 1 in a field twice higher. 
Nevertheless, even such a small Stark shift is still within 
the range of energies that are detectable with current 
experimental techniques (see, e.g. Ref. [3^. It is also 
clear that electric fields applied in x-direction affect more 
the exciton polarization and, consequently, lead to higher 
Stark shift as compared to a field applied in ^-direction. 

Another interesting feature of the exciton Stark effect 
in black phosphorus, as compared to the one in conven¬ 
tional bulk semiconductors, is the opportunity to observe 
such effect also in an excited state. The same study we 
made for the ground state exciton in Fig. [2] was also done 
for the second excited state {n = 3), where quite simi¬ 
lar features are observed, although small deviations from 
quadratic Stark shift are observed at high applied fields, 
which also occurs for carbon nanotubes exciton Stark 
shift. @ A clear difference, however, lies in the stronger 
Stark shift observed in the n = 3 case, as compared to the 
previous case, which is due to its lower binding energy: 
in the = 4 case, for instance, a AE meV shift is 
obtained for an electric field of 50 kV/cm applied in the 
x-direction. A similar shift is also obtained for n/ = 1, 
but for a higher field, = 110 kV/cm. On the other 
hand, such lower binding energy also forbids us to inves¬ 
tigate the effect up to higher fields without dissociating 
the electron-hole pair. The polarizabilities of the exci¬ 
ton ground and second excited states, for systems with 
different number of layers, are summarized in Table nn 
where the anisotropy of the Stark effect becomes quite 
evident by the distinct polarizability parameters found 
for electric fields applied in different directions. 

The polarization induced by the field is expected to 
contribute to a broadening of the excitonic peak in the 
absorption spectrum and, therefore, hinder the visualiza¬ 
tion of the Stark effect. Nevertheless, our results for the 
oscillator strength in Fig. [3] suggest that such contribu¬ 
tion is quite small in the ground state, especially for a 
monolayer: |'0i(O,O)p changes by only 2% (1%) for 
fields up to iAfy) = 210 kV/cm, whereas such a change 
is already observed for a field of Fx = 90 {Fy = 70) in the 
case of four layers. As for the excited state, |'03(O,O)p 
changes up to ^ 14(27)% in the case of Fy = 50 (30) 
kV/cm for ni = 1(4), as compared to the zero field case. 
Results in each panel are divided by the zero field os¬ 
cillator strength of the ni = 1 case. This means that 
the ni = 2,3 and 4 systems exhibit ground (second ex¬ 


cited) state oscillator strengths that are ^ 55.5% (61.9%), 
38.5% (46.7%), and 28.8% (36.3%) of that of the mono- 
layer case, respectively. Notice that the influence of the 
electric field on the ground state oscillator strength is 
higher when the field is applied in the x-direction (closed 
symbols), whereas for the n = 3 state, it is higher for F 
applied in the ^-direction (open symbols). This can also 
be understood as a manifestation of the anisotropy of the 
exciton wave function: although n = 1 and 3 are both 
s-like states, we observe that (x^) > {y‘^) for the former, 
whereas < {y‘^) for the latter (cf. supplementary 
material), which explains such opposite behaviour of the 
oscillator strengths of these states under applied fields. 
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FIG. 3: Oscillator strengths for the (a) ground (n = 1) and 
(b) second excited (n = 3) exciton states as a function of 
the strength F of an electric held applied in the x (closed 
symbols) and y (open symbols) directions. The symbols refer 
to the different number of layers, as in Fig. 1. 


In conclusion, we have calculated the excitonic states 
in few layer black phosphorus in the presence of an ex¬ 
ternal electric held, by numerical diagonalization of the 
effective mass Hamiltonian. The method developed here 
is shown to be easily adapted for other layered mate¬ 
rials and agrees well with recent experiments on WS 2 . 
Due to the non-Coulomb effective interaction potential 
between electron and hole in such 2D system, the de¬ 
generacy coming from the SO(3) symmetry of the pla¬ 
nar hydrogen atom model is lifted, whereas the bands 
anisotropy is responsible for lifting the degeneracy of an¬ 
gular momentum eigenstates. Even so, the exciton spec¬ 
trum in black phosphorus can still be satisfactorily htted 
by a hydrogen-like expression. Due to the high binding 
energies, excitons can withstand strong in-plane electric 
helds without dissociating, This allowed us to ob¬ 
serve Stark shifts up to ^15 meV, with electric fields up 
to ^200 kV/cm, not only for the ground state exciton, 
but also for excited states, without significant deprecia¬ 
tion of the oscillator strength. We believe that the clear 
and unusually anisotropic exciton Stark shift predicted 
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here may stimulate further photoluminescence and re¬ 
flectance experiments in few layer black phosphorus un¬ 
der in-plane electric fields in the near future. 
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In order to find the eigenstates of the Hamiltonian Eq. 
(3) of the main manuscript, we propose a finite differ¬ 
ence scheme on a variable mesh, which allows us to take 
better account of the singularity of the electron-hole in¬ 
teraction potential as the relative coordinate r ^ 0. The 
x^-plane is discretized in x Ny points (here, we take 
Nx = Ny = 1020) separated by = x{y)j — x{y)j-i^ 
where the mesh varies according to -h 

{1 - exp [-(j - A^3.(y)/2 - l)/(5®(s^)]}. This sets 

the values of for j > N^^y)/2, namely, for x{y) > 0. 
The grid is then mirror reflected to the negative part of 
the x(^)-axis. In this way, we can control how close the 
points are to r = 0 by adjusting tune how nar¬ 

row the mesh is in the neighborhood of this point by 
adjusting and set a larger mesh far from r = 0 by 

adjusting Results in this paper are obtained with 

YYu = 5 X 10-5 A, = 100 and = 0.41. 

After discretization, we write the two-dimensional 
Hamiltonian as a five-diagonals matrix and numerically 
diagonalize it. The grid points, labeled in the (line, col¬ 
umn) form (i,j), are then mapped into a single label 
1. For a non-uniform mesh, diagonalizing this matrix is 
cumbersome, because the discretized Hamiltonian matrix 
is asymmetric: 


Hu = 


1 1 1 It. 

+ ~rXTo -^ Tv - To -^ 1 V tO -^ 
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Hi,1+1 = - 


Hi 


l,lTNy 
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Hl,l-Ny = — 
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where = {h^jf)\i + ) /2- Nevertheless, 

symmetrization of the Hamiltonian can be achieved by 
multiplying all the terms in the Hamiltonian by •, 

or equivalently, B = jCjCH, where £ is a diagonal 
matrix whose diagonal elements are LxjLy^i. Since 
jCjCHjC~^ jCf; = ECCfj, then = ECCfj or, 

equivalently C~^BC~^^ = E^ (where ^ = Cijj). Hence, 
diagonalizing C~^BC~^, which is symmetric, gives the 
same eigenenergies as diagonalizing the Hamiltonian it¬ 
self. This matrix is then numerically diagonalized by 
Arnoldi’s method in order to obtain the exciton states. 

The exciton energies obtained from this method are 
shown in Fig. 1 of the main manuscript. As for the 
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FIG. 1: (a) Average width of the exciton function in each 
direction for m — 1 (black circles) and 4 (green stars). Open 
(full) symbols are results for ((x^)). (b) Wave functions 
for the four low-lying exciton states, considering = 4 black 
phosphorus in the absence of electric fields. 


eigenfunctions. Fig. [T]^a) shows the width of each exci¬ 
ton state wave function, considering ni = 1 (black cir¬ 
cles) and 4 (green stars). Linear interpolating curves are 
shown just as a guide to the eye. The anisotropy (namely, 
the difference between (open symbols) and {y‘^) (full 
symbols)) is more evident for some specific exciton states. 
As verified by the wave functions for the four low-lying 
exciton energy states in n/ = 4 case, shown in Fig. m, 
the n = 1 state resembles a Is orbital, as expected, but 
its width is twice larger in the direction of lighter effec¬ 
tive mass, namely, the x-direction. Subsequent states 
resemble 2py, 2s and 2px orbitals, respectively. Similar 
wave functions are also found in Fig. 7 of Ref. [l| for 
ni = 1. The n = 2 state exhibits ^ (^^), which 
is unusual for a Py state, that is normally more spread 
in the ^-direction. However, this is compensated in the 
case of phosphorene by its lighter effective mass in the 
x-direction. The n = 3 state, which is a 2s-hke orbital. 
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exhibits a slightly higher width in the ^-direction, which 
is more evident in the ni = A case, as compared to n/ = 1. 
This has consequences for its oscillator strength, as dis¬ 
cussed in the main manuscript. From the widths shown 
in Fig. HJa), one can infer that the exciton wave function 
in both directions is spread over tens of angstrons, i.e. 
over a region many times larger than the inter-atomic 
distance a ^ 2k in black phosphorus, Q as required for 
the validity of the Wannier-Mott theory of excitons used 
here. In fact, just like in carbon nanotubes, Hi exci¬ 
tons in few layer black phosphorus undergo low screen¬ 
ing outside the layer, resulting in quite large electron-hole 
binding energies, as opposed to those of Frenkel excitons. 
Nevertheless, the dielectric function in the layer itself is 
on the same order of magnitude of those of other semi¬ 
conductor materials, [H such as Si and GaAs (namely, 
e ~ lOeo), which leads to the large spatial extension of 
the exciton wave function observed here, thus justifying 


the use of Wannier-Mott theory. 
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